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1.   INTRODUCTION 

The  recent  development  of  several  digital  computer  programs  for 
the  analysis  of  electrical  circuits  has  added  a  new  dimension  to  pro- 
cedures for  designing  networks.  These  network  analysis  programs  have 
the  capability  to  calculate  time  and  frequency  response  data  from  a 
topological  description  of  an  electric  circuit.   The  response  data 
generated  can  be  used  by  the  designer  to  alter  the  network  parameters 
and  configuration  until  a  satisfactory  design  is  achieved.   If  this 
approach  to  network  design  is  followed,  it  is  important  to  have  a 
digital  computer  operating  mode  with  minimal  turn-around  time  so  that 
the  designer  can  observe  the  effects  of  parameters  changes  and  proceed 
to  an  acceptable  design  in  an  efficient  manner. 

An  alternative  approach  is  to  use  a  digital  computer  program  to 
make  design  decisions.   When  this  is  done,  a  network  analysis  program 
supplies  performance  data  which  is  used  to  alter  the  network  parameters 
in  a  systematic  way.   The  parameters  are  adjusted  to  minimize  a  per- 
formance measure;  the  performance  measure  (or  objective  function)  is 
a  function  which  indicates  the  deviation  of  the  network  response  from 
the  desired  response.   Minimization  of  the  performance  measure  (or 
optimization)  can  be  carried  out  by  any  of  several  computational 
algorithms. 

Regardless  of  whether  the  design  is  performed  by  repeated  analysis, 
or  by  optimization,  the  computer-aided  approach  offers  several 
advantages  over  classical  design  procedures,  e.g.  Darlington  synthesis, 
Brune  cycle,  etc.    In  computer-aided  design  procedures,  parameter 
constraints,  nonlinearities,  and  parasitic  effects  can  be  included. 
The  design  of  a  network  of  specified  complexity  whose  response  approxi- 
mates a  desired  response  can  also  be  carried  out.   In  addition,  compu- 
tations can  be  performed  using  the  physical  design  parameters,  such  as 
transistor  base  width,  or  a  film  thickness  in  an  integrated  circuit,  as 


the  variables.   References  1-4  provide  a  good  summary  of  the 
state  of  the  art  of  computer-aided  network  design. 

The  investigation  reported  here  has  been  limited  to  the  design 
of  linear  networks  in  the  frequency  domain.   Three  optimization 
algorithms  have  been  used  in  combination  with  the  network  analysis 
program  CALAHAN  to  solve  several  example  problems, 


2.   NETWORK  DESIGN  AS  A  NONLINEAR  PROGRAMMING  PROBLEM 


The  basis  for  computer-aided  network  design  by  the  use  of 
numerical  optimization  procedures  is  the  formulation  of  the  design  as 
a  problem  in  nonlinear  (or  mathematical)  programming.   A  nonlinear 
programming  problem  is  one  in  which  a  function  J  (the  objective 
function)  of  r  variables  p  ,  p^,  ...,  p   is  to  be  minimized  (or 
maximized)  subject  to  constraining  relations  of  the  form 


g^P^  •••>  Pr)  >  0   ,  i  =  l,2,...,m   . 


(1) 


To  simplify  the  notation,  vector-matrix  symbols  will  henceforth  be 


used,  that  is,  p 


_ ) 


,  and  (1)  can  be  written 


g.(p)  >  0  ,  i 


-  1,2, ... ,m  , 


(la) 


or,  simply 

g(p)   ^  0  , 


(lb) 


where 


g(p)  - 


gx(p) 
g2(P> 


The  set  of  points  which  satisfy  all  of  the  constraints  given  by  (1)  is 
called  the  feasible  set.   The  constraints  in  (1)  that  are  satisfied 
with  the  equality  holding  are  said  to  be  active. 
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The  feasible  set  is  assumed  to  be  a  closed,  bounded,  and  convex 
set,  and  the  objective  function  is  assumed  to  have  continuous  second 
partial  derivatives  with  respect  to  the  components  of  p. 

In  the  network  design  application,  the  objective  function  J 

provides  a  measure  of  the  deviation  of  the  network  frequency  response 

from  the  desired  frequency  response,  and  the  variables  are  a  set  of  r 

network  parameters  p  ,  p„,  ...,  p  .   In  general,  the  performance 

measure  will  be  assumed  to  have  the  form 
N 
J(p)  =   E  h(G(f.),  GD(f.),  f.)  (2) 

i=l 

where  f.(Hz),  i  =  1,2,...,N,  is  a  set  of  frequencies  for  which  the  desired 
response  is  specified,  G  (f.)  is  the  desired  response  at  f  =  f., 
and  G(f.)  is  the  response  of  the  network  at  the  frequency  f..   The 
scalar  function  h  satisfies  the  following  properties: 

1.  h  >  0,  for  all  G(f  ),  G  (f  )  ,  f    ,   i  =  1,2,... ,N     (3a) 

2.  h  =  0,  only  if  G(f.)  =  Gff.)  (3b) 

l     D  l 

In  the  example  problems  discussed  subsequently,  the  performance 
measured  used  was 

N 
J(p)  =  E   w.[20  4og10  |GD(f.)|  -  20  4og1()  |G(f.)1]       (4) 
i=l 

where  the  w.  are  positive  numbers  called  weighting  factors.   By 
inspection  of  (4)  it  is  observed  that  the  larger  a  particular  w.  is, 
the  larger  the  value  of  J  for  a  given  response  deviation;  hence,  by 
adjusting  the  w.'s  the  actual  response  can  be  made  to  conform  more 
closely  to  the  desired  response  at  certain  frequencies  at  the  cost  of 
larger  deviations  at  other  frequencies. 

In  this  investigation  the  constraints  considered  were  restricted 
primarily  to  the  form 

3i  "  Pi  ~   bi  '   i  =  l>2>"->r      •  (5) 


where  a   and  b.  represent  the  lower  and  upper  bounds  for  the  ith 
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parameter.   These  constraints  can  be  put  into  the  form  of  (1)  by 
defining 


§1^> 

— 

PL    '    31 

>   0 

g2(p) 

= 

"Pl  +  bl   " 

■   0 

82r-l(P>  =   Pr  -  ar  *  ° 

g2r(p)    =   "Pr  +  br  >  0  (6) 

For  some  of  the  optimization  algorithms  considered,  extension  to  other 
types  of  constraints  is  a  routine  matter;  this  will  become  evident  in 
the  subsequent  discussion. 


3.   SOLUTION  OF  NETWORK  DESIGN  PROBLEMS 

To  design  a  network  by  using  the  procedures  described  here  the 
designer  must  first  select  the  configuration  of  the  network.   Some 
investigators  have  reported  the  incorporation  of  logic  into  design 
programs  for  "growing"  network  elements,  but  this  investigation  did  not 
include  this  possibility. 

After  specifying  the  network  configuration,  a  performance  measure 
and  any  weighting  factors  must  be  selected.   The  desired  frequency 
response  characteristic  must  be  specified  either  in  the  form  of  a  table 
of  values,  or  as  an  analytical  expression  that  can  be  included  in  the 
program. 

Once  these  preliminary  steps  in  the  design  procedure  have  been 
completed,  and  the  necessary  data  has  been  provided  for  the  combined 
analysis-optimization  program,  the  design  parameters  are  adjusted  by 
the  computer  program  until  the  performance  measure  is  minimized.   The 
network  analysis  program  provides  frequency  response  data  to  the 
optimization  program  which  then  adjusts  the  variable  parameters  in  a 
systematic  way  until  the  minimum  is  found.   The  nature  of  the  data 
that  must  be  provided  by  the  analysis  program  depends  on  the  optimiza- 
tion algorithm  being  used.   Generally,  the  network  analysis  program 
will  be  required  to  provide  either  values  of  J  only,  or  values  of  J 
and  the  first  partial  derivatives  of  J  with  respect  to  the  elements  of 

A  flow  chart  of  the  computational  procedure  is  shown  in  Figure  1. 
It  should  be  noted  that  although  this  flow  chart  and  the  procedure 
outlined  are  for  frequency-domain  design,  they  apply  aiso  to  design 
in  the  time  domain  by  replacing  the  work  "frequency"  by  the  word  "time", 
The  time  response,  however,  would  be  calculated  within  the  network 
analysis  program  in  a  different  manner  than  the  frequency  response. 


Yes 


Initialize  the  programs  and  read 
input  data  including  the  initial 
guess  for  the  variable  parameters 


Adjust  the  network  parameters 


Use  the  network  analysis  program 
to  calculate  the  frequency  response 
of  the  network 


Print  out  optimum  parameter 
values,  graph  response  if  desired 


Figure  1   Flowchart  of  the  computational  procedure 


4.   THE  NETWORK  ANALYSIS  PROGRAM 

There  are  many  network  analysis  programs  that  have  been  written. 
Most  of  these  programs  are  available  from  the  authors  or  from  the 
sponsoring  agencies.   A  few  of  the  more  widely  known  analysis  programs 
are  CALAHAN,  SCEPTRE,  ECAP,  CORNAP,  and  NASAP.   In  this  investigation 
the  program  CALAHAN  has  been  used  exclusively  because: 

1.  It  is  well  documented  and  hence  easy  to  use. 

2.  It  provides  both  time  and  frequency  response  data  for  linear 
networks  (including  networks  with  dependent  sources). 

3.  The  program  is  compatible  with  the  IBM  360/67  in  use  at  the 
Naval  Postgraduate  School. 

Experience  has  indicated  the  CALAHAN  is  adequate  for  use  with  the 
optimization  programs  employed;  however,  it  is  clear  that  modifications 
could  be  made  which  would  make  the  program  operate  more  efficiently 
when  used  in  an  iterative  mode.   This  comment  will  be  amplified  after 
reviewing  the  use  of  the  CALAHAN  network  analysis  program. 

To  use  CALAHAN,  a  topological  description  of  the  network  and 
numerical  values  for  its  parameters  must  be  provided  as  input  data  to 
the  program.   The  network  to  be  analyzed  is  considered  to  be  in  the 
form  of  a  two-port  black  box  as  shown  in  Figure  2.   One  of  the  user 
options  in  CALAHAN  is  the  selection  of  the  network  function  of  interest; 


+ 


vL(t) 


o 


Linear 
Network 


i2(t) 
I 0  + 


v2(t) 


Figure  2   Two-port  network 


the  choice  is  made  by  specifying  the  input  parameter  KEY  1.   The  alter' 
natives  are  shown  in  Table  1. 


KEY  1 

Network  Function 

Symbolic 

1 

Voltage  transfer  function 

V2 
Vl 

V-o 

2 

Open  circuit  driving  point 
impedance 

r2-° 

3 

Open  circuit  transfer 
impedance 

V2 

ia-o 

4 

Short  circuit  driving  point 
impedance 

Xl 
Vl 

v2-o 

5 

Short  circuit  transfer 
admittance 

Vl 

v2  =  0 

6 

Current  transfer  function 

v2  =  o 

TABLE  1.   Network  function  options 

The  other  input  data  that  must  be  provided  to  the  modified  version  of 
CALAHAN  is: 

1.  The  number  of  passive  (RLC)  elements. 

2.  The  number  of  controlled  sources  (current  sources  that  are 
voltage  controlled). 

3.  The  number  of  nodes. 

4.  The  positive  input  node. 

5.  The  negative  input  node. 

6.  The  positive  output  node. 

7.  The  negative  output  node. 

8.  KEY  1 


9.  The  location,  element  type,  and  value  of  each  of  the  network 
parameters . 

10.  The  frequency  response  data. 

Items  1  through  9  are  placed  on  a  single  data  card  according  to  the 
format  (12,  IX).   Item  9  is  provided  by  supplying  one  data  card  for 
each  element;  this  data  card  contains  the  node  numbers  at  the  ends  of 
the  element,  the  element  type,  and  the  element  value  according  to 
the  format  (2(12,  IX),  Al,  IX,  F10.0).   The  information  for  Item  10  is 
provided  on  two  cards:  the  first  card,  which  is  read  according  to  the 
format  (II,  IX,  13) ,  contains  a  1  or  2  in  column  1  depending  on  whether 
a  linear  or  logarithmic  frequency  scale  is  desired,  and  columns  3 
through  5  contain  either  the  total  number  of  frequency  intervals 
(for  a  linear  scale),  or  the  number  of  points  per  decade  (for  a  log 
frequency  scale);  the  second  card,  which  is  read  according  to  the 
format  (2F10.0),  contains  the  smallest  and  largest  frequencies. 
Figure  3  illustrates  a  circuit  diagram  and  the  input  data  to  compute 
the  frequency  response  V2(f)/V.(f)  for  f  =  0.150  Hz.  to  f  =  0.250  at 
intervals  of  \f  =    .001  Hz.   The  "output"  data  that  result  in  this 
example  are  the  values  of  V  (f)/V  (f)  for  f  =  0.150,  0.151,  ...,  0.250. 
The  unmodified  version  of  CALAHAN  provides  additional  data,  including 
the  poles  and  zeros  of  the  network  function,  and  plots  of  magnitude 
and  phase  versus  frequency;  however,  since  this  information  is  not 
desired  when  the  analysis  program  is  used  in  an  iterative  mode,  it 
has  been  eliminated  from  the  program. 

As  mentioned  previously,  the  use  of  CALAHAN  with  an  iterative 
optimization  algorithm  is  somewhat  inefficient.   The  reason  for  this 
is  that  each  time  CALAHAN  is  called,  the  tree-generating  subroutine 
begins  as  if  a  new  network  were  being  analyzed.   As  a  result,  a 
disproportionate  amount  of  time  is  spent  generating  information  that 
could  be  computed  once  and  stored. 
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rrr\ 


V. 


"3 

_fY-r\ 


5 


R 


v. 


L   =  1.5451  H. 
L„  =  1.3820  H. 


Lr  =  0.3090  H. 


C2  =  1.6944  F. 
C.  =  0.8944  F, 
R   =  1.0  ft 


(a)  The  circuit  diagram 


|  Columns  |  Columns  |  Columns 
!    1  -  10  I   11  -  20  I  21  -  30 


Card  #  9 
Card  #  8 
Card  #  7 
Card  #  6 
Card  #  5 
Card  #  4 
Card  #  3 
Card  #  2 
Card  #  1 


J 1 1 

fb.150             b.250 
1                        1 

/  11    100                                       | 

1 

^[05    02    R  1.0 

! 

/  04    05    L   0.13090 

|                           | 

/j04    02   C   0.18944 

'    03    04    L    l.|3820                 j 

'  !03    02   C    lj.6944 

1                          1                           ' 

^!01    03    L    lj5451                 ' 

^06    00    05    dl    02    05    02:  01 

(b)  The  input  data 
Figure  3 
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5.   THE  OPTIMIZATION  ALGORITHMS 

Three  optimization  algorithms  have  been  combined  with  CALAHAN 
(in  separate  programs)  to  perform  optimal  design  of  linear  networks 
in  the  frequency  domain.   A  description  of  these  optimization 
algorithms  follows. 

5 . 1   The  Gradient  Projection  Method 

The  gradient  projection  algorithm,  developed  by  J.  B.  Rosen 
(see  reference  5)  determines  the  minimum  of  a  function  of  several 
variables  which  are  constrained  to  lie  in  a  closed  and  bounded  convex 
region  R  in  the  parameter  space  P.   R  is  defined  as  the  interior  of  a 
convex  region  bounded  by  a  set  of  linear  constraints.   A  detailed 
exposition  of  Rosen's  method  is  given  in  reference  5  and  also  in 
reference  6.   Here,  a  brief  description  of  the  important  features  of  the 
algorithm  will  be  given. 

The  algorithm  is,  as  the  name  implies,  a  gradient  search,  but 

with  the  capability  of  solving  problems  in  which  the  variables  are 

constrained.   The  version  of  Rosen's  algorithm  that  has  been  used  in 

this  investigation  assumes  that  the  constraints  are  linear,  that  is, 

T 
N,P  -  v  £  0   ,  (7) 

where  N.  is  a  known  r  x  I   matrix,  p  is  the  vector  of  r  variable 
parameters,  v.  is  a  known  vector  of  I   constants,  and  the  superscript  T 
denotes  matrix  transposition.   By  selecting  I   =  2r  and 


•5*  = 


&  i  4] 


an 

a2 

Zi  = 

a 

r 

-bi 

"b2 

• 
• 

-b 
r 
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the  formulation  in  (7)  yields  the  special  case  given  in  (6). 

The  performance  measure  J  is  minimized  subject  to  the  constraints 
given  by  (7).   The  strategy  employed  in  the  search  procedure  is  to  move 
in  the  direction  of  steepest  descent  (the  negative  gradient  direction) 
in  the  parameter  space  so  long  as  the  function  values  decrease  and 
none  of  the  constraints  are  violated.   If  a  move  in  the  negative 
gradient  direction  would  cause  one  or  more  of  the  constraints  to  be 
violated,  the  negative  gradient  is  projected  onto  the  constraint 
boundary  at  the  point  where  the  constraint  violation  would  occur. 
This  projection  of  the  negative  gradient  is  accomplished  by  a  q  x  q 
projection  matrix  P  ;  q  is  the  number  of  constraints  in  (7)  that 
are  active  (satisfied  with  the  equality). 

The  projection  matrix  is  formed  from  the  columns  of  the  matrix 
N  which  correspond  to  active  constraints.   The  active  constraints  can 
be  written  as 

where  v   is  formed  by  selecting  those  entries  of  v.  which  correspond 
to  active  constraints.   The  projection  matrix  is  given  by 

P„  =  I  -  N  [N  T  N  ]_1  N  T  (9) 

~q    „        ~q  ^q  ^q-J    _q 

The  algorithm  is  carried  out  in  such  a  manner  that  P  can  be  computed 

by  using  recurrence  equations  that  do  not  require  inversion  of  the 

T 
matrix  TN   N  1; this  feature  provides  a  tremendous  computational 

~q  ~qJ 

saving.   If  the  current  best  point  does  not  lie  on  the  boundary,  but 

within  R,  then  P  is  the  identity  matrix. 

(i)    ^ 
If  p    is  the  current  best  point  the  next  parameter  vector  is 

given  by 

p(i+D  =  p(i)  +  ^q[-oJ(X(i))/^p]/||SJ(^p(i))%||  (10) 

where  j    (the  step  size)  is  calculated  by  using  an  interpolation 
procedure  or  a  single  variable  search. 
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To  illustrate  how  the  gradient  projection  algorithm  operates,  a 
simple  example  is  shown  in  Figure  4.   In  this  example  the  constraints 
are  given  by 


o  <  ?i   <  1 


0  *  p2  ^  1   , 


(ID 


or 


1  0 


0   10-1 


(11a) 


Several  equal-value  contours  of  the  function  J  are  shown  in  Figure  4 
and  R  is  the  region  bounded  by  these  constraints.   Noffeice  that  R  is  a 
convex  set  because  every  point  on  a  straight  line  drawn  between  any 
two  points  in  R  is  also  in  R. 


If  p' 


is  the  first  trial  point,  then  z     is  the  unit  vector  in 


the  negative  gradient  direction  and  a  move  is  made  to  the  boundary  point 
p^  '    on  the  line  H, .   At  p^   ,  moving  in  the  negative  gradient  direction 
would  violate  the  constraints,  so  a  move  is  made  in  the  direction  of 
the  vector  z     (a  unit  vector  in  the  direction  of  the  projection  of 
the  negative  gradient  at  p    onto  the  line  H, ) .  Moving  in  this 
direction  yields  p    where   further  motion  in  the  direction  of  z 
would  violate  the  constraint  p_  <  1.   Repeating  the  steps  used  earlier 

provides  the  information  that  no  improvement  can  be  made  by  moving  from 

(2) 
pv  '\    hence,  the  iterative  procedure  terminates. 

This  algorithm  has  been  coded  in  FORTRAN  IV  by  the  principal 
investigator;  in  reference  7  a  description  is  given  of  the  combined 
operation  of  CALAHAN  and  the  gradient  projection  algorithm. 

Two  features  of  the  CALAHAN  -  gradient  projection  program  are  of 
special  interest. 

1.   The  gradient  projection  algorithm  requires  not  only  values  of 
J,  but  also  the  partial  derivatives  of  J  from  which  the 
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J(i)  <  J<k>.  for  i  >k 


Figure  4   Geometric  interpretation  of  the  gradient  projection 
algorithm 
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negative  gradient  vector  is  formed.   These  partial  derivatives 

were  approximated  by  a  finite  difference  approximation,  that 

is 

SJ        J(P1,P2,. • • ,Pj+A, ...,Pr)  -  J(P1,P2» • • • »Pj» • • • >Pr) 

5^  (£>  " I 

Using  this  approach,  (r+1)  network  analyses  are  required  to 
evaluate  the  gradient  vector  at  the  point  p. 
2.   The  interpolation  scheme  for  finding  the  step  size  t  suggested 
in  reference  5  requires  additional  evaluations  of  the  gradient 
vector.   Because  the  gradient  calculation  is  costly  in  terms 
of  the  number  of  calculations,  an  alternative  approach  was 
adopted  in  which  only  the  values  of  J  were  needed.   In  addition 
to  proving  more  efficient  in  this  application,  this  single 
variable  search  procedure  alleviated  difficulties  believed  to 
be  caused  by  the  occurrence  of  local  minima. 
Reference  7,  in  addition  to  describing  the  combined  operation  of 
the  programs,  also  provides  several  examples  which  illustrate  appli- 
cations . 

5.2   The  Pattern  Search  Method 

Unlike  the  gradient  projection  method,  which  requires  the  first 
partial  derivatives  of  J  with  respect  to  each  of  the  variable  parameters, 
the  pattern  search  algorithm,  reference  8,  requires  only  values  of  the 
function  J.   Perturbations  in  the  elements  of  the  parameter  vector  p 
are  used  to  determine  the  direction  of  a  change  in  p  which  decreases 
the  value  of  J.   The  perturbation  data,  which  provides  local  information, 
is  then  used  to  suggest  larger  changes  in  p  which  should  further 
decrease  J.   These  larger  changes  in  p  are  called  pattern  moves.   The 
algorithm  alternates  perturbation  moves  (or  local  explorations)  with 
pattern  moves.   Each  time  a  successful  pattern  move  occurs,  the  step 
size  of  the  following  pattern  move  is  increased;  thus,  continued 
successful  pattern  moves  cause  future  moves  to  become  larger  and  larger. 
When  a  pattern  move  eventually  fails,  the  pattern  is  destroyed  and  the 
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algorithm  attempts  to  begin  generating  a  new  pattern.   The  search 
terminates  when  an  exploration  move  cannot  decrease  the  function,  and 
the  perturbation  step  size  has  been  reduced  to  a  value  below  a 
pre-selected  minimum. 

Figure  5  illustrates  a  simple  two-parameter  problem  whose  solution 
indicates  the  steps  followed  by  the  pattern  search  algorithm.   The 
figure  shows  the  two-dimensional  parameter  space  and  several  equal- 
value  contours  of  J. 

At  the  initial  trial  point,  b    ,  the  function  J  is  evaluated,  and 
this  value  is  denoted  by  J(b    ).   The  parameter  p   is  perturbed  by  an 
amount  +  Ap   and  the  function  value  increases,  indicating  a  failure, 

so  p1  is  perturbed  by  -  APi  and  the  function  value  decreases  (a  success) 

(1) 
from  its  value  at  b    .   Perturbation  of  p  by  an  amount  Ap9  further 

~  (2) 

decreases  the  function  and  yields  b    as  the  current  best  point. 

Using  the  rationale  that  further  perturbations  would  likely  result  in 

the  same  pattern  of  successes,  the  next  move  is  to  the  point  t    , 

where 

t<X>  =  b(1>  +  2[b<2)  -  b(1)]  .  (13) 

Notice  that  the  function  is  not  evaluated  at  t    unless  all 
perturbations  about  t    produce  values  of  J  greater  than  J(b   )  -- 
which  does  not  occur  for  the  geometry  of  this  example.   The  iterative 
procedure  continues  in  this  manner  with  perturbation  cycles  and  pattern 
moves  alternating.   The  pattern  moves  are  determined  by  using  the 
equation 

_t(i)  =b(1)  +2[b(1+1)  -b(1)]  (14) 

Several  additional  moves  are  also  shown  in  Figure  5. 

As  mentioned  previously,  the  pattern  moves  continue  to  increase 
in  length  as  long  as  successes  occur.   If  a  failure  occurs  on  a  pattern 
move  [all  perturbations  about  t    produce  values  of  J  larger  than 
^(l+l)^  and  j(t(L))  >  j(b^1+  '))],  then  the  pattern  is  destroyed. 
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J(i>  <  J(k> 


for  i  >  k 


(2) 


Figure  5   Geometric  interpretation  of  the  pattern  search  algorithm 
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When  it  becomes  impossible  to  make  successful  perturbation  moves  from 
a  current  best  point,  the  perturbation  step  size  is  decreased  by  a 
specified  amount.   When  the  perturbation  step  size  has  decreased  below 
a  pre-specif ied  level,  the  iterative  procedure  terminates. 

The  pattern  search  method  is  available  as  subroutine  DIRECT  in  the 
Naval  Postgraduate  School  Computer  Facility  FORTRAN  Library.   In 
reference  9  a  description  is  given  of  the  combined  operation  of  CALAHAN 
and  the  pattern  search  algorithm.   The  version  of  pattern  search  reported 
in  reference  9  was  modified  to  incorporate  constraints  of  the  form 

a.  <  p.  ^  b.    ,   i  =  1,2,... ,r  (15) 

in  the  algorithm.   Reference  9  also  presents  several  network  design 
examples . 

5.3   The  Fletcher-Powell  Method 

The  Fletcher-Powell  method,  reference  10,  is  basically  a  gradient 
method;  however,  unlike  other  gradient  techniques,  the  convergence  of 
the  algorithm  becomes  more  rapid  as  the  minimum  is  approached.   The 
reason  for  this  is  that  the  gradient  vectors  calculated  in  the  iterative 
procedure  are  used  to  generate  an  estimate  of  the  matrix  of  second 
partial  derivatives.   As  the  procedure  evolves,  the  estimated  matrix 
of  the  second  partial  derivatives  approaches  the  actual  matrix  of 
second  partial  derivatives  and  the  convergence  of  the  algorithm  becomes 
quadratic. 


The  motivation  for  this  approach  is  provided  by  considering  the 
>r  series  expansi 
of  up  to  second  order 


Taylor  series  expansion  of  J  about  the  point  p    which  includes  terms 


-  P(i)] 


<E(i+1))  ■  <.(1))  ♦  [| G(1))J    i 

+  \ [p(1+1)  -  z(i)]  T  fi  da))  [i(i+1)  "  £(l)' 


(16) 
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where 


hJ 


bJ 


dPj    BP2 


*P, 


denotes  the  gradient  of  J  with  respect  to  p,   —  (p   )  is  the 

2 
gradient  vector  at  the  point  p    and  — s-  (p    )  is  the  matrix  of 

second  partial  derivatives 


*Pi' 


s2j 


b2j 


^PrSpL 


32J 


BPLBP2 


Bp2^P1     Sp2' 


a2j 


BPrdp2 


(i) 


,2J 


BpL^Pr 

82J 
9p2sPr 


B2J 


*P, 


evaluated  at  the  point  p 

The  point  p    is  known  (initially  it  is  a  guess)  and  the  point 
p      is  to  be  found.   The  gradient  of  J  with  respect  to  p      is 


dj/  (i+l)> 
3FU    J 


*P 


BJ  f    (i)\    S  J  MO 


+ 


^P 


P         "  P 


,   (17) 


,  ..   (i+1)  .  .  .  .  .      ... 

and  if  p      is  to  be  the  minimizing  point,  it  is  necessary  that 


I  ^(i+0)  = 


=  0  . 


(18) 
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Hence, 


P        "  P 


1  fAzW)'^W 


(i+1) 


Q"1  grad  J(i) 


(19) 


This  equation  suggests  that  to  find  p      the  change  in  p  should  not 
be  in  the  negative  gradient  direction  (this  would  be  the  case  if  the 
matrix  of  second  partials  were  the  identity  matrix),  rather  a  deflected 
gradient  should  be  used.   If  the  matrix  of  second  partials  Q  (the 
Hessian  matrix)  were  known,  Equation  (19)  could  be  used  as  the  basis 
for  an  algorithm  having  quadratic  convergence.   Unfortunately,  when 
an  explicit  relationship  for  J  in  terms  of  p  is  not  known,  as  will  be 

the  case  in  computer-aided  network  design  applications,  the 

2     2 
computation  of  ^  J/dp   is  very  difficult,  and  attempting  to  find  the 

matrix  of  second  partial  derivatives  by  perturbations  often  leads  to 

inaccurate  results. 

In  the  Fletcher-Powell  Method,  the  gradients  themselves 
are  used  to  obtain  an  estimate  of  Q.   For  a  function  J  that  is  quadratic 
in  p  (equation  (16)  is  exact),  the  algorithm  converges  in  r  iterations, 
where  r  is  the  dimension  of  p. 

The  iteration  equation  used  is 


p(i+D    =   p<i>-   T   tt<i>    grad   j<« 


(20) 


where  t  is  the  step  size  which  is  determined  by  performing  a  single 
variable  search.   The  single  variable  search  in  the  subroutine  being 
used  involves  first  an  extrapolative  phase,  during  which  the  value  of 
t  which  makes  J      the  smallest  along  the  line  in  the  deflected 
gradient  direction  is  bracketed,  and  an  interpolation  phase,  during 
which  the  minimizing  value  for  t  is  determined. 

The  Fletcher-Powell  Method  as  described  is  capable  of  solving 
unconstrained  minimization  problems;  however,  in  the  class  of  problems 
of  interest  here  the  variable  parameters  will  be  required  to  satisfy 
constraints.   To  allow  application  of  the  Fletcher-Powell  Method  some 
way  of  incorporating  the  constraints  must  be  found.   One  way  of 
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handling  constraints  is  to  use  penalty  functions;  the  performance 
measure  J  is  altered  by  some  additional  terms  which  depend  on  the 
constraints.   Intuitively  speaking,  the  terms  in  the  altered  performance 
measure,  J  ,  that  correspond  to  the  constraints  should  be  large  if  the 

3 

constraint  boundaries  are  being  approached  and  small  if  the  parameters 
are  within  the  region  where  the  constraints  are  satisfied. 

To  be  specific,  assume  that  the  constraints  are  of  the  form 

g.(p)  >  0,   i  =  1,2,. ..,m.  (21) 

As  usual,  it  is  desired  to  minimize  J.   The  method  used  is  to  form  the 
altered  performance  measure  (devised  by  Fiacco  and  McCormick,  reference 

11) 

1 
J  (P)  -  J(P)  +  R  S   — 7-r  .  (22) 

s~    ~    i=L  gi(.p; 

R  is  a  positive  constant. 

If  the  initial  point  lies  in  the  interior  of  the  feasible  region, 
i.e.,  where 

g.(p)  >  0,   i  =  1,2,. ...td   ,  (23) 

then  if  the  iterative  procedure  adjusts  p  to  a  value  near  the 
boundary,  J   approaches  infinity  because  at  least  one  of  the  g.  terms 
approaches  zero.   In  this  manner  the  penalty  term 

1 
R   E 


1-1   gi(P} 

tends  to  keep  the  values  of  p  within  the  admissible  region.   If  R  is 
allowed  to  become  very  small,  and  the  constrained  minimum  lies  on  the 
boundary,  the  value  of  p  found  by  minimizing  J  will  approach  the  true 
minimum  point  p*. 

To  illustrate  the  Fiacco-McCormick  procedure,  the  following 
example  is  presented.   The  function  to  be  minimized  is 
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J(P)  -  P  (24) 

and  the  constraints  are 

1  <  P  ^  2  .  (25) 

Putting  this  constraint  in  the  form  of  equation  (21)  gives 


gL(p)  =  p-1  >  0 
g2(P)  =  2-P  '■-   °  » 


hence 


J   =  p  +  R  , 
a         I  p-1 


+ 


2-P 


(26) 


(27) 


Plots  of  J   as  a  function  of  p  for  two  values  of  R  are  shown  in 
a  r 

Figure  6. 

Generally,  the  minimization  procedure  is  performed  several 
times  with  decreasing  values  of  R.   It  can  be  shown  (reference  10) 
that  doing  this  generates  a  sequence  of  points  which  converges  to  p* 
(provided  that  J  meets  certain  mathematical  requirements). 

3. 

In  combining  the  Fletcher-Powell  and  Fiacco-McCormick  procedures 
with  CALAHAN  it  has  been  assumed  that  the  constraints  are  of  the  form 
given  in  equation  (6),  thus, 


and 


J  (p,R)  =  J(p)  +  R  S 

j=l 


hJ 


1 


p. -a. 
J   J 


+ 


1 


•p.+b. 
J   J 


(28) 


OP. 


(i) 


+  R 


-1 


1 


(P 


(i) 


V 


c-p^  +  V2 


(29) 


where  p.  is  the  jth  component  of  the  vector  p. 

The  Fiacco-McCormick  approach  has  been  combined  with  the  Naval 
Postgraduate  School  Computer  Facility's  version  of  the  Fletcher- 
Powell  algorithm  (subroutine  FLAP)  and  CALAHAN.   Several  modifications 
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Figure  6   The  augmented 
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have  been  made  to  improve  the  operation  of  the  combined  program.   One 
of  these  modifications  ensures  that  p  is  never  adjusted  to  lie  out- 
side of  the  feasible  region.   Such  an  adjustment  can  occur  in  the 
Fletcher-Powell  method  during  the  extrapolation  phase  of  the  procedure 
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6.   ILLUSTRATIVE  EXAMPLES 

In  this  section  two  numerical  examples  will  be  presented.   The 
three  optimization  techniques  discussed  previously  were  applied  to 
these  examples  to  provide  a  basis  for  comparing  the  algorithms. 
Although  these  examples  have  networks  with  only  passive  elements,  the 
programs  are  also  capable  of  optimizing  linear  networks  containing 
active  elements.   References  7  and  9  present  several  additional 
example  problems,  including  some  witi  non-ideal  network  elements. 

Example  1.   The  problem  is  to  determine  the  reactive  element  values 
in  the  network  shown  in  Figure  7  so  that  the  frequency  response 


+  1  o- 


L_nnr\ 

L0 


4  _rm 

Lr 


-o  + 


2o- 


1.0  Q 


Figure  7.   Fifth-order  low  pass  filter 


minimizes  the  performance  measure 

10 
J  =   S  [20  jtog10|GD(f.)|  -  20  £og10|G(f.)l]    . 
i  =  l 


(30) 


G(f)  =  Vc-^  (f )  /V.  ~  (f )  is  the  actual  frequency  response  and  the  desired 
frequency  response  characteristic   is  given  by 


w  = 


vTiT 


(31) 


(2TTf.) 


2N 


where  N  =  5.   This  desired  response  is  the  Butterworth  response 
characteristic  for  a  fifth-order  filter;  if  the  element  values 
in  the  network  shown  in  Figure  7  are 
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L  =  1.545  H 


1.694  F    L^  =  1.3820  H   C,  =   0.8944  F 
3  4 


Lj.  =  0.3090  H   the  desired  response  characteristic  will  be 
obtained  exactly.   The  frequencies  used  were  f   =  .15  Hz.,  .16  Hz., 
...  ,  .24  Hz. 

The  results  obtained  are  summarized  in  Table  2.   Although  the 


Gradient 

Pattern 

Search 

Projection 

Fletcher 

- Powe 1 1 

Fifth-Order 

Initial 

Final 

Initial 

Final 

Initial 

Final 

Butterworth 

Value 

Value 

Value 

Value 

Value 

Value 

Value 

Ll 

1.545 

1.600 

1.540 

1.000 

1.557 

1.600 

1.475 

C2 

1.694 

2.000 

1.700 

2.400 

1.478 

2.000 

1.800 

L3 

1.382 

1.500 

1.380 

1.100 

1.695 

1.500 

1.312 

C4 

0.894 

0.900 

0.895 

0.900 

1.064 

0.900 

0.877 

L5 

0.309 

0.500 

0.307 

1.000 

0.500 

0.500 

0.242 

Table  2   Element  values  for  fifth-order  low  pass  filter 

element  values  obtained  by  the  pattern  search  and  Fletcher-Powell 
algorithms  differ  significantly  from  the  Butterworth  values,  the 
frequency  response  curves  obtained  are  almost  identical  with  the  desired 
frequency  response  curve.   The  frequency  response  results  obtained  from 
the  gradient  projection  method  are  shown  in  Figure  8. 

The  Fletcher-Powell  method  took  eighteen  iteration  cycles  to 
determine  the  values  given  in  Table  2.   [An  iteration  cycle  is  defined 
as  a  minimization  of  J  for  a  specified  value  of  the  factor  R  in  Equa- 
tion (22)].   If  the  algorithm  operated  as  anticipated,  one  iteration 
cycle  should  have  been  sufficient  to  obtain  convergence.   The 
slowness  of  convergence  is  attributed  to  the  relative  insensitivity 
of  the  performance  measure  to  variations  in  the  parameter  values; 
this  means  that  near  the  minimum  the  gradients  calculated  by  the 
perturbation  method  are  likely  to  be  too  inaccurate  to  obtain  the 
expected  convergence  properties  of  the  algorithm. 
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Figure  8   Actual  and  desired  frequency  responses 
for  low  pass  filter 


Example  2.   The  reactive  element  values  in  the  network  shown  in  Figure  9 
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Figure  9   Bandpass  filter  configuration 
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are  to  be  determined  to  minimize  the  performance  measure 

J  =   S   [20  iog10|GD(f.)|  -  20  £og10|G(f.)|]' 

where  G(f)  -  V2(f)/V1(f)  and  20  Xog1Q|G  (f .) |  is  the  bandpass  filter 
response  shown  in  Figure  10;  a  tabulation  of  the  gain  values  on  the 
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Figure  10   Actual  and  desired  bandpass 
filter  frequency  responses 


desired  frequency  response  curve  is  contained  in  reference  9. 

The  element  values  obtained  by  the  three  algorithms  are  contained 
in  Table  3. 

The  performance  of  Fletcher-Powell  method  was  less  successful 
than  in  Example  1.   The  values  in  Table  3  were  obtained  in  one 
iteration  cycle,  but  correspond  to  a  performance  measure  of  J  =  175.8. 
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Gradient 

Pattern 

Search 

Projection 

Fletcher 

-Powell 

Initial 

Final 

Initial 

Final 

Initial 

Final 

Value 

Value 

Value 

Value 

Value 

Value 

Ll 

0.500 

0.278 

1.000 

0.231 

0.500 

0.193 

cl 

1.200 

0.968 

1.200 

1.140 

1.200 

1.213 

L2 

0.750 

0.512 

1.000 

0.520 

0.500 

0.433 

C2 

0.750 

0.500 

1.000 

0.485 

1.200 

0.988 

L3 

0.500 

0.231 

1.000 

0.237 

0.750 

0.050 

C3 

1.200 

1.070 

1.200 

1.140 

0.750 

0.737 

Table  3   Bandpass  filter  element  values 

By  contrast,  the  performance  measure  which  corresponds  to  the  pattern 

-fci 
search  values  is  of  the  order  of  1  x  10 

The  frequency  response  curve  obtained  for  the  gradient  projection 

element  values  is  shown  in  Figure  10. 
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7.   CONCLUSIONS 

Three  optimization  techniques  have  been  combined  with  the 
network  analysis  program  CALAHAN  and  used  to  solve  several  illustrative 
examples.   Although  the  examples  contained  only  networks  with  passive 
elements,  the  programs  are  capable  of  optimizing  linear  networks  with 
active  elements.   Of  the  three  methods  pattern  search  seems  to  be 
most  effective;  however,  it  is  believed  that  this  conclusion  is 
mainly  due  to  the  finite  difference  approximation  used  for  the 
gradient  in  the  other  two  algorithms.   The  finite  difference 
approximation  is  also  believed  to  be  responsible  for  the  difficulties 
encountered  with  the  Fletcher-Powell  algorithm. 

To  improve  the  efficiency  of  the  computational  procedures  several 
modifications  could  be  made.   First,  the  network  analysis  program 
should  be  altered  so  that  once  the  appropriate  network  trees  are 
calculated,  they  should  be  stored  instead  of  being  re-computed  each 
time  the  analysis  program  is  called;  this  change  should  reduce  the 
required  computation  times  significantly.   Second,  the  gradient 
vector  should  be  calculated  by  using  formulas  rather  than  by  using 
a  finite  difference  approximation.   This  modification  should  make 
the  gradient  projection  and  Fletcher-Powell  algorithms  far  more 
efficient.   Finally,  investigation  of  techniques  for  altering  the 
network  configuration  as  well  as  its  element  values  should  broaden 
the  class  of  problems  for  which  these  computer-aided  design  methods 
are  useful. 
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